NONEQUILIBRIUM SHEAR VISCOSITY COMPUTATIONS WITH 

LANGEVIN DYNAMICS 
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Abstract. We study the mathematical properties of a nonequilibrium Langevin dynamics which 
can be used to estimate the shear viscosity of a system. More precisely, we prove a linear response 
result which allows to relate averages over the nonequilibrium stationary state of the system to equi- 
librium canonical expectations. We then write a local conservation law for the average longitudinal 
velocity of the fluid, and show how, under some closure approximation, the viscosity can be extracted 
from this profile. We finally characterize the asymptotic behavior of the velocity profile, in the limit 
where either the transverse or the longitudinal friction go to infinity. Some numerical illustrations 
of the theoretical results are also presented. 
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1. Introduction. The determination of the macroscopic properties of a mate- 
rial given its microscopic description is the fundamental goal of statistical physics [1] . 
Macroscopic properties can be classified into two categories: (i) equilibrium proper- 
ties, such as the heat capacity or the equation of state of the system (relating the 
pressure, the density and the temperature); and (ii) transport properties, such as the 
thermal conductivity or the shear viscosity. The determination of transport prop- 
erties is conceptually and numerically more challenging than the determination of 
equilibrium properties since transport phenomena depend both on the chosen ther- 
modynamic ensemble and on the prescribed microscopic dynamics (which has to leave 
the thermodynamic state of the system invariant). Often, the Hamiltonian dynamics 
is considered as the reference dynamical evolution of the microscopic system. How- 
ever, this dynamics exactly preserves the energy of the system, while energy exchanges 
with the environment are expected to happen. We believe that the choice of the un- 
derlying dynamics of the system is a modelling choice. In any case, a careful study of 
the dependence of the computed transport properties as a function of the parameters 
of the dynamics should be performed. 

There are roughly three types of methods for computing transport coefficients: 

(i) Transient methods, where the relaxation of some local disturbance is monitored 
as a function of time, and macroscopic coefficients are obtained by fitting the 
observed response to the evolution predicted by a macroscopic evolution equa- 
tion. For instance, a local temperature hot spot can be created in the middle 
of a homogeneous material, and, assuming that the heat equation describes well 
the evolution of the kinetic temperature field, the diffusion of the energy allows 
to estimate the thermal conductivity of the material (see for instance [El [31] ) ; 

(ii) Equilibrium methods, which are based on time integrals of correlation functions, 
the so-called Green-Kubo formulas [8j [20] . These correlation functions are ob- 
tained by sampling initial conditions according to the thermodynamic ensemble 
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at hand, and averaging over all possible evolutions from these initial conditions. 
For example, the shear viscosity of a fluid is given by 



where \T>\ is the volume of the simulation box, o~ xy (t) is the off-diagonal term of 
the microscopic stress tensor at time t, T is the temperature of the system, and 
(•) denotes a canonical average; 
(iii) Nonequilibrium methods, where the system is in a steady-state characterized 
by the existence of stationary fluxes and spatial gradients of some quantities. 
These fluxes and gradients can be controlled by forcing terms acting on the 
boundaries of the system (for instance, a velocity profile can be obtained by fixing 
the average velocity in the extremal slabs of a fluid, see the early review [15], 
or subsequent works such as [5]), or by a bulk process where fictitious forces 
act on all particles (the so-called synthetic molecular dynamics approach [B]). 
Boundary driven methods are usually numerically less efficient than bulk-driven 
methods since there are more correlations in the system and the convergence to 
a steady-state is slower. 
A review of the most standard approaches to compute the shear viscosity can be read 
in [6]. See also |33j for a focus on nonequilibrium methods. 

We focus in this study on bulk-driven nonequilibrium molecular dynamics tech- 
niques. Since the system is driven out-of-equilibrium by a non-gradient force, some 
thermostatting mechanism is required to prevent the uncontrolled increase of the 
energy and to ensure that a steady-state can indeed be reached. In many works fo- 
cusing on the computation of shear viscosity, the thermostatting is performed with 
deterministic dynamics, such as Nose-Hoover like thermostats [MJ [2] or isokinetic 
dynamics (see [6j Chapter 3]). The ergodicity of these dynamics is at most unclear 
(and non-ergodicity can even be proved rigorously in some limiting cases |21l 122]). 
The mathematical analysis of these methods is therefore untractable. Only formal 
results of linear response can be written down. In some studies, the thermostatting 
is performed using dissipative particle dynamics (301 119j , which includes stochastic 
terms. The ergodicity of these dynamics is however a very difficult issue, and the only 
existing results we are aware of concern one-dimensional systems |29) . 

We decided to use a standard Langevin dynamics as the underlying dynamics of 
the system since this dynamics is ergodic and has many nice mathematical properties, 
while still being close enough to the Hamiltonian dynamics. In essence, the nonequi- 
librium dynamics we propose is obtained from the standard Langevin dynamics by 
adding a nongradient force, which can be interpreted as some fictitious external forc- 
ing term. The effect of this term is to create a velocity profile in the direction of the 
forcing, and the viscosity of the fluid can be extracted from this profile. The novelty 
of this work with respect to the (numerous) existing studies on the computation of 
shear viscosity is the rigor of the mathematical arguments used to prove linear re- 
sponse results and obtain the effective equation on the observed velocity profile in 
terms of the applied external force. In particular, we benefited from recent develop- 
ments on hypocoercivity [36]. Besides, one of our main concern is the dependence of 
the viscosity as a function of the parameters of the underlying dynamics, in particular 
the friction. We analyzed the large friction asymptotics by extending and adapting 
mathematical studies of the auto-diffusion coefficient jTOJ [571 E] • F° r the low friction 
dependence, we rely on numerical simulations. 
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This paper is organized as follows. We start by describing the Langevin dynamics 
we propose in Section[3J and show the existence and uniqueness of the stationary state. 
The mathematical properties of this stationary state are studied in Section|3](with the 
proofs postponed to Section [5]). In particular, we give a rigorous proof of the linear 
response. We then show how to compute the viscosity, and characterize its asymptotic 
behavior for large frictions by determining the limiting behavior of the velocity profile. 
We finally present some numerical illustrations of the theoretical results in Section 21 

2. The nonequilibrium Langevin dynamics. 

2.1. Description of the dynamics. We consider a system of iV particles, en- 
closed in a periodic simulation box T>. For simplicity of notation, we restrict ourselves 
to two-dimensional systems, so that V = L X T x L y T (T being the one-dimensional 
torus R/Z). The particles are described by their positions q = (qi, . . . , (jfjv) € V N 
and their momenta p — (pi, . . . ,pn) G R 2A \ and we assume that they have identical 
masses m > 0. Our results can however straightforwardly be extended to the gen- 
eral case of particles with different masses and/or three-dimensional systems. We 
write q, = (q xl ,q yi ) G £>, pi = {p xi ,p y i) G M 2 , as well as q x = (q xl , . . . , q xN ), 
Px = (Pxi: ■ ■ ■ ,Pxn) and similar definitions for q y ,p y - The mass density of the system 
is 

mN 

\T>\ = L x L y being the volume of the box. Finally, we denote by 

n 2 

»=i 

the Hamiltonian of the system, the function V being the potential energy. 

The equations of motion we propose are a linear perturbation of the Langevin 
equations, with some additional non-gradient external force in the x-direction (the 
direction of the flow). The dynamics reads (for i = 1, . . . , N): 

aqt.t = — at, 
m 



dp xi<t = -V fei Vfe) dt + iF{q y%t ) dt - lx ^ dt + J^f dW t\ (2.1) 



dp y%t = -y qyi V{q t ) dt - ly ^ dt+J-& dW t 



m 



where ft — l/(fcsT) is the inverse temperature, ^ is the magnitude of the external 
nongradient force, (WjF, Wf )t>o is a 2iV-dimensional standard Brownian motion and 
the friction coefficients "f Xl "f y are real positive numbers. In order to avoid irrelevant 
technical issues, we make the following assumption: 

Assumption 1. The potential V and the external force F belong respectively to 
C°°(V N ) and C°°(L y T). 

Note that the canonical measure 

M<l,P)dqdp = Z- 1 C -P H ^ dqdp (2.2) 

is invariant by the dynamics (|2.1[) when £ = 0, and is actually the only invariant 
measure (see for instance the discussion in [3]). 
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2.2. Existence and uniqueness of an invariant measure. When £ ^ 0, 
there is no obvious invariant probability measure, and the very existence of such a 
measure is not guaranteed a priori. However, in the case when 7^,7^ > 0, standard 
techniques based on Lyapunov functions and hypoellipticity arguments can be resorted 
to to prove the existence and uniqueness of an invariant measure which has a smooth 
density with respect to the Lebesgue measure. For further purposes, it is convenient to 
consider the reference space L 2 (ipo) (where the measure tpo is defined (|2.2p ). endowed 
with the scalar product 



(f>9) L »(ii> a ) := / f{v,p)g{q,p)ipo{q,p)dqdp. 

JV N xR 2N 

We can then state the following result. 

Theorem 1. Consider 7 x ,7a > and suppose that Assumption^ holds. Then, 
for any the dynamics (|2.1[) has a unique smooth invariant measure with density 

i>{. e C°°(V N x R 2N ). Besides, there exists £* > such that, for any £ G 
the following expansion holds in L 2 (ipo): 

1>t = fvh, /e = i + E^' ( 2 - 3 ) 

where f& £ L 2 (ipo) is such that ||ffc||i, a (i/io) 5~ C(£*) _fc . 

The proof is presented in Section al I The existence and uniqueness of an invariant 
measure in the case when either 7^ = or 7^ = is a much more difficult question. 
To obtain such a result, more precise assumptions on the potential are required (see 
Remark [1] in Section POT) . 

In the sequel, averages with respect to the measure ip£ are denoted 



We = / h{q,p)ijj i (q,p)dqdp= (h,^)^^). 

JV N xR 2N 

3. Mathematical analysis of the viscosity. 

3.1. Linear response. Linear response results allow to compute the average of 
some properties with respect to the nonequilibrium measure in terms of equilibrium 
averages, in the limit when the parameter giving the strength of the nonequilibrium 
forcing vanishes. To describe the result more precisely, we introduce the infinitesimal 
generator associated to the equilibrium Langevin process (i.e. (|2.ip in the case when 

e = o)= 



-^0 — ^ham "^Uhm; 



where 



and 



Aam = — • Vg - VV{q) • V p 

m 



Ahm = 5^ 7a [~ ■ V Pa + 4 A Pc «) = ^j- Tadiv Pa (e^V^-) . 

a — x,y ^ ^ a—x^y 

It can be proved (see Section UTTj) that Aq 1 is a well defined operator on the Hilbert 
space 



n = \fe L 2 (^ ) 



M> = } = n {i}- 

v n xr 2N 
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where the orthogonality is with respect to the -L 2 (i/>o) scalar product. 
Define also the adjoint operator on L 2 (ipo) of the generator: 

The generator of the nonequilibrium perturbation reads 

N 



i=i 

and its adjoint on L 2 (ipo) is 

N 8 
^— ' m 

i=l 

The generator of the dynamics (|2.ip is therefore 

Az = Ao + iB, 

with adjoint L^ = Lo + £,B* . 

Linear response is an easy consequence of Theorem [TJ 

COROLLARY 1. Under the same assumptions as in Theorem [7J and for any 
function h regular enough, 

(Aoh) { 8 I ^ \ 



1=1 I £ 2 «>o) 



Besides, for any function h £ W, 

N 



5 \ »=1 / L^ ) 



Proof. Since (A h, l)m^„) = (h, L 1} L 2 (i/)q) = 0, it holds 

(A h) i = (Aoh,h)m^) = £MoMi>£»(*o + °(£ 2 ) = Mi W o) + o(e 2 ). 

Now, in the proof of Theorem [I] we show that (see (|5.ip ) 

L f! = -B*l = --J2p xi F(q yi ), 

which gives the expected result. □ 

3.2. Local conservation of the longitudinal velocity. We prove in this sec- 
tion a conservation equation for velocities in the ^-direction, when spatial averages 
over small windows in the transverse direction y are considered. This allows to state 
an equation relating the off-diagonal term of the stress tensor and the nongradient 
force acting on the system, see (|3.6p below. Our derivation may be seen as a mathe- 
matically rigorous counterpart to the seminal work of Irving and Kirkwood [17] . We 
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assume from now on that the potential energy is given by a sum of pairwise interac- 
tions: 

V( qi ,...,q N )= Yl V (lft-^l), (3-2) 

l<i<j<N 

for some given smooth potential V. 

Consider the following average longitudinal velocity: 

T N 

L£(Y, q,p) = Y,P*iXe (q yi - Y) , (3.3) 
1=1 

where Xe (with < e < 1) is an approximation of the identity on L y T. More precisely, 

/ % 1 is — nL,, 



where \ G C°°(R) has support in [0, L y ] and J V X = 1- The factor L y 
accounts for the fact that \s has units of inverse lengths: in fact, 



is the average velocity of the system. In practice, averages such as (|3.3[) are computed 
with bin indicator functions (see Section |4.1[) . 

We also need a spatially localized (with respect to the altitude Y) version of 
the off-diagonal term of the stress tensor. This quantity is given by the following 
expression (the fact that it can be interpreted as some stress tensor is motivated 
below by the limiting spatial average (|3.5|) as well as the conservation law (|3.6|) ): 



x V i=l \<i<j<N w Hn J 1v3 

Note that the limiting spatial average over Y 



lim-1 f \ xy (Y,q,p)dY 



N 



1 | V~~^ PxiPyi V~~^ mUI |\ Wxi Qxj)(qyi Qyj) 



(3.4) 



(3.5) 



L X L V I *-r i m ' & — o, 

• ' \ ; i i<i<j<N HJI 

is the standard expression encountered for the off-diagonal term of the pressure tensor 
without spatial localization. The expression (|3.4j) comes out naturally from the math- 
ematical analysis (see the proof of Proposition [lj , and was already proposed in [34] 
(where it is called the 'method of planes'). 

The relationship between the local longitudinal velocity and the off-diagonal term 
of the stress tensor is made precise in the following proposition. 

Proposition 1. The limits 

u r (Y) = lim lim - 

v ; e->0£-vO £ 
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and 



belong to C°°(L y T) and 



zyK 1 ! = 11111 11111 7 



o~xy(Y) = lim lim 



^p. +7x p Ux (Y)=pF(y), (3-6) 

where ~p — p/m is the particle density. 

The proof is based on an application of Corollary[T]with fl3.3[) as a test function h. 
The order of the limits e — > and £ — > cannot be inverted since the linear response 
result of Corollary [T] cannot be applied with h replaced by the limit of \e (which is a 
Dirac mass). 

Equations similar to (|3.6p could be written down for other quantities such as the 
transverse velocity U y and longitudinal and transverse energy fluxes (see [17] for the 
original derivation of the corresponding equations). 

3.3. Definition and closure relation for shear viscosity computations. 

We now discuss a closure relation for (|3.6[) . which allows to obtain an equation on the 
average velocity only, from which the viscosity can be extracted. 

By analogy with continuum fluid mechanics, we define the shear viscosity rj as 
follows: 

a xy (Y):=- V (Y)^±. (3.7) 
This definition leads to the following equation on u x : 

-4? ( V(Y) ^p-) + l- P u x {Y) = pF(Y). 



dY V dY 

In bulk homogeneous fluids, the simplest closure is to assume that 

r)(Y) =r)>Q, (3.8) 



so that the following equation on u x is obtained: 

- nu'^Y) + lx -pu x {Y) = pF(Y). (3.9) 

In order to ensure the uniqueness of the solution when j x = 0, an additional condition 
should be added (such as a vanishing integral over the domain L y T). 

The equation Q3.9P obtained with the help of the closure relation is the basis for 
numerical methods to compute the shear viscosity given a potential energy function V. 
We were not able to justify mathematically the assumption Q3.8p . We nonetheless 
provide a numerical validation of this assumption in Section [4.2.21 

3.4. Asymptotic behaviour of the viscosity for large frictions. An im- 
portant issue is the dependence of the viscosity on the parameters of the dynamics. 
For the Langevin dynamics (|2.1[) . this means understanding the dependence of the 
viscosity on the friction parameters 7 X , j y . The limits 7 X — > or ^ y — > are very diffi- 
cult to study mathematically without strong assumptions on the potential and/or the 
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geometry of the system (see Remark [1} . We therefore rely on numerical simulations 
for these cases (see Sections 14.2.31 and I4.2.4[) . 

On the other hand, the limit when one of the friction parameters goes to infinity 
can be studied. To this end, we have to understand the limit of the velocity field u x 
as either j x or j y goes to infinity. This is done by rigorous asymptotic analysis. 
Thanks to ()3.9|) , limiting behaviors of the viscosity may be inferred from the limiting 
behaviors of the velocity profiles. The key result to obtain the limiting velocity profile 
is to characterize the limit of some averages with respect to specific solutions of the 
Poisson equation (see (|3.11|) and p.!3[) below). 

3.4.1. Infinite transverse friction. We start with the case 7 y — > +oo, for a 
fixed value j x > 0. 

Theorem 2 (Infinite transverse friction) . Consider a given smooth function G 
and a longitudinal friction j x > 0. Define Ao^y) := .4o = -4ham + j x A x . thm + 

JyAy,thm, With 

A a>t hm = -— ■ V Pq + \ A Pa , (3.10) 

m p 

and denote by / 7 the unique solution in % of the equation 

N 



Then, there exist f , f 1 G ^(ipo) an d a constant C > such that, for all j y > j x , 

ll^-/ -^ 1 / 1 !^^^^ (3.12) 
Besides, the function f is of the general form 

N 

f°(q,P) = ^2G{q yl )<p l (q x ,q y ,p x ), 

i=l 

where the functions <pi are C°° . 

The proof can be read in Section \5. 31 The above result can be used to understand 
the limit of u x (Y) as 7 a — > +oo. Indeed, by Proposition [U 



u^(Y) := lim { 3^)k = _A (j2p xi F(q yi ),%l(Y,q,p) 



\i=l 



where —Ao{^ y )^ y {Y,-) = U X (Y,-) is a Poisson equation of the form (|3.11[) (with 
G(y) proportional to Xs(y — Y)). The convergence result (13. 12)) shows that ^/^(Y, •) 
has a limit as j y — > +oo, and the limiting velocity field reads 

<'%Y) = |^ (^P x iF(q yi ),J2xe(q yj -Y) ( f> j (q x ,q y ,p x )\ 

The latter quantity has a limit as e —> 0, so that the velocity field converges to 
some limiting field Therefore, the viscosity extracted from (|3.9I) also has a finite 
limit. These theoretical considerations are illustrated by numerical simulations in 
Section 
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3.4.2. Infinite longitudinal friction. We now consider the limit ~/ x — > +00, 
for a fixed value j y > 0. In this case, the leading term of the expansion in inverse 
powers of j x is 0, and a refined convergence result is needed to discuss the limit of 
the velocity profile. 

Theorem 3 (Infinite longitudinal friction) . Consider a given smooth function G 
and a transverse friction ~/ y > 0. Define Aq{^ x ) ■— Ao = Aham+JxA x ,thm+"/ y A y ,thm, 
and denote by / 7x the unique solution in H of the equation 

N 

- A(7=e)/ 7i , = y^ J PxiG{q y i). (3.13) 

Then, there exist f 1 , f 2 G H 1 (iPq) and a constant C > such that, for all 7^ > j yi 

C 

ll 

Besides, the dependence of the function f 1 in the variable p x can be written explicitly 



N 



l^-^f-^fWHHM^-i- ( 3 - 14 ) 



f l {q,p) = m^2p xl G(q yl ) + f (q,p y ). 



i=i 

The proof can be read in Section [5^1 To obtain asymptotics on the velocity field, 
we apply the above convergence result with G(y) proportional to Xe{y — Y) (denoting 
by fl the first term in the expansion in inverse powers of 7^): 

S ' \i=l / £»(V>o) V ' x 

PL I N N \ ( 1 

/ 22p xi F(q yi ),22p xj Xe(q v j ~ Y) ) + 0' 



Nmy.. 



J- 1 ' L>ty 



^ [ " F(y) Xe (y-Y)dy + o(^ 

ix Jo \ii 



where we have used the fact that {p X i, fl) L 2 (tp ) = ^ smce fl does not depend on p x . 
This shows that the following limit is well defined: 

u x {Y) = \im lim lx u%{Y) = F(Y). (3.15) 

The limiting velocity profile u x does not depend on the specific interaction potential V, 
and is the same for all systems with pairwise interactions. Besides, the viscosity n 
cannot be extracted from Q3.9P since (u £ x )" is of order 7" 1 while F and ^i e x u x are of 
order 1. The limit 7^ — ► +00 is therefore somewhat degenerate from a theoretical 
viewpoint. Numerical simulations however allow to investigate the large 7 X asymp- 
totics, see Section 14.2.41 

4. Numerical results for the Lennard-Jones fluid. We present in this sec- 
tion some numerical illustrations of the theoretical results obtained in Section [3J 

4.1. Numerical implementation. 
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Description of the system. We consider a Lennard- Jones fluid, which is a standard 
test case for shear flow computations, in a 2-dimensional setting (in order to limit the 
number of degrees of freedom and henceforth obtain results with lower statistical 
uncertainties) . The potential energy is of the form (|3.2p , with 




V LJ (r)=4, LJ N^j -^fj j. (4.1) 

Actually, it is numerically more convenient to work with a truncated potential, which 
reads: 

if T ^spline? 

V( r ) = { v spiinc (r) if r splinc < r < r cut , 




if r > r, 



cut ■ 



The function w sp iino is a polynomial of order 3 which is such that the potential is C 1 
on (0, +oo) (note that there is a singularity at r = so that this potential does not 
satisfy Assumption [1] However, it seems that this singularity does not show up in 
the numerical simulations. Any problem related to this singularity could be overcome 
by modifying appropriately the potential for the very small values of r). We use 
r C ut = 3d L j and r sp i inc = 0.9r cut . 

All the results presented below are in reduced units, which are determined by 
setting to 1 the energy £lj, the length Jlj, and the mass m. The remaining tunable 
parameters of the model are the force amplitude £ and the friction parameters j x , j y . 

The thermodynamic state of the system is determined by the fluid mass density p 
and the temperature T. In the numerical illustrations below, we set /J = 0.4, p — 0.69 
and consider L x — 360 and L y = 18. The number of simulated particles is therefore 
N = 4500. 

We have checked that the thermodynamic limit is attained for the systems we 
simulate, i.e. that the values of the viscosity and the profiles we present are converged 
with respect to increasing values of L x ,L y (at fixed density), see [T5] . 

Nongradient forces. We consider three different external perturbations, which are 
all normalized so that — 1 < F(y) < 1: 

/ 2iry\ 

(i) sinusoidal perturbation: F(y) = sin I I • 

(ii) piecewise linear perturbation: F(y) = 



(iii) piecewise constant constant perturbation: F{y) = 




-I. • .'/ /.,. 



Note that only the sinusoidal force satisfies Assumption[TJ This shape of perturbation, 
introduced in [7], is the most popular choice for shear viscosity computations. 

Integration of the dynamics. The dynamics (|2.1|) is discretized using a standard 
splitting scheme, similar to the schemes proposed in |23j . The evolution is decomposed 
as the superposition of (i) a Hamiltonian part, which is integrated with the standard 
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Verlet scheme [351 15]: and (ii) a fluctuation/dissipation part containing also the non- 
gradient force, which can be integrated analytically since it is an Ornstein-Uhlenbeck 
process with a constant drift. The numerical scheme reads 

At 



q n+1 = q n + Atp n 



-1/4 



1+1/2 _ p n+l/4 _ ^yy^n+l) 



(4.2) 

where a XtV = exp(— 7 X)V Ai), and G",G™ are independent and identically distributed 
standard Gaussian random variables. Note that this scheme is well behaved in the 
limits j x — > +oo and/or ^ y — > +oo, as well as in the limits "/ y — > or j x — ¥ (the 
well posedness of the latter case is a consequence of the limit (1 — a x ) /j x — » At as 
lx -+ 0). It reduces to the standard Verlet scheme when F = and r ) x —ly = 0. 

We use At = 0.005 in all the simulations below. This time step ensures that the 
relative error in energy is about 1% for the Verlet scheme. 

Numerical localization. To analyze the various fields which can be constructed 
from the numerical data generated by the simulation (longitudinal velocity, off-diagonal 
component of the stress tensor, kinetic temperatures, etc), we use a binning procedure 
in the Y variable. More precisely, we introduce a mesh with a uniform spacing AY, 
centered on the altitudes Y s = (s + 1/2) AY (with < s < S - 1 and SAY = L y ). 

The microscopic observables we wish to average are either the longitudinal ve- 
locity (|3.3p or the off diagonal stress tensor Q3.4p . Both functions are of the general 
form 

N 



A s (Y,q,p) = Y,a i (q,p)$ s (q vi -Y), 



i=i 

where <!> e is either Xe or an integral of this function. Averages of such functions over 
each cell are computed as 

Taking advantage of the integration in the Y variable, it is possible to take the limit 
e — > in the latter expression. This amounts to compute ensemble averages with 
respect to bin indicator functions (or their integrals). For instance, the average lon- 
gitudinal velocity in the sth bin is 

^ = £NmAY \ ^ Pxil [Ys-AY/%Y s +AY/2](qyi)) ■ 

In practice, the ensemble average (-)j is computed as a time average over trajectories 

(« n ,P n )»=l,..,JW 



12 



R. Joubaud and G. Stoltz 



Estimation of the viscosity. The solutions of Equation (|3.9[) are periodic in the 
y-variable and are hence most easily analyzed using Fourier series (see for instance 
the discussion in [51 H2]). We consider the Fourier coefficients of the the average 
longitudinal velocity u x and the force F, given respectively for k S Z by 



1 [ L » . . /2ifc7ry\ , 1 [ Ly . . / 2ifc7ry 

u k = Y~ u x (y)exp ' jay, F k = — I F(y) exp 

^y Jo V / ^j/ JO 



^y 



dy. (4.3) 



The coefficients Uk can be esimtated numerically using trajectory averages as 



- a : i , ■ 1 V s - ™~ f 2lk7Tq yj 



r^LL-^Uf • («) 

n— 1 j — 1 v y / 

This is a valid estimation provided the marginal distribution in the position of one 
particle is the uniform law on the domain. By translation invariance, this is true when 
no external force is present. It remains approximately true when £ is not too large. We 
checked that this approximation has no influence on the presented numerical results. 

The shear viscosity is obtained from a Fourier analysis of (|3.9|) . The value of r\ 
should be independent of k € Z. It should also satisfy the following equation: 

U k = f * . (4.5) 



P \ L y 

The shear viscosity is finally obtained as 



U k 'J \2kw 

A confidence interval on the value of Uk can straightforwardly be obtained from 
the estimator (|4.4[) using block averaging procedures. The statistical uncertainty on 
the viscosity is then obtained with (14.61) . Since the coefficients Uk decrease very 
rapidly as |fc| increases, the relative statistical errors increase rapidly as well. We 
therefore restricted ourselves to \k\ = 1 in our numerical simulations. 

4.2. Numerical results. In all cases, time averages were computed over iVjter — 
10 7 iterations. 

4.2.1. Linear response. We first verify numerically the linearity of the lon- 
gitudinal velocity as a function of the magnitude £ of the nongradient force. More 
precisely, we check that \U\\ is constant, for the three forces F at hand, in the case 
when (j x ,Jy) — (1,1) (see Figure |4~T|) . In the sequel, unless otherwise stated, the 
numerical results are obtained with £ = 0.1. 

The numerical results show that linear response is valid even for values of £ large 
enough. In fact, a more refined analysis (see |18j ) shows that, even if no nonlinear 
effect can be observed on the longitudinal velocity for the values of £ we considered, 
nonlinear effects on the kinetic temperature cannot be ignored for values of |£| > 0.1. 

Other numerical simulations (not reported here, see \18\ ) show that linear response 
actually applies in the case when 7 X = although this case is not covered by our 
theoretical analysis. 
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Fig. 4.1. Value of \U\ as a function of £, for ("fxi'fy) = (1,1) and the three nongradient forces 
at hand (piecewise constant perturbation X; piecewise linear perturbation +; sinusoidal perturbation 
o)- 

4.2.2. Validation of the closure. We present in Figures T4.21 14.31 and 14.41 the 

numerical approximations of the longitudinal velocity u x and the off-diagonal term of 
the stress tensor a xy . The latter function is compared to the quantity —T)u' x , where rj 
is obtained from (|4.6[) . and u x is evaluated using a second order finite difference. The 
good agreement between a xy and —r\u' x validates the assumption (|3.8p . the discrepan- 
cies resulting from statistical fluctuations magnified by the numerical derivative, and 
also, for the piecewise constant force, from the singularity at L y /2. 

Besides, the velocity profile is consistent with (|3.9I) (as can be checked by com- 
paring the numerical solution and the solution of ()3.9j) computed with the value of 7/ 
estimated from the simulation). 




value value 



Fig. 4.2. Velocity profile and off diagonal component of the stress tensor for the sinusoidal 
nongradient force. 

4.2.3. Dependence and asymptotics in j y . We first verify numerically that 
the velocity profiles converge to some limiting profile, and in fact that U\ converges 
to some limiting value U£° (see Figure |4~5|) . We estimated by long simulations 
with a y = in (14. 2[) (which amounts amounts to formally setting ~f y to +oo), and 





Fig. 4.4. Velocity profile and off diagonal component of the stress tensor for the piecewise 
constant force. 



computed the distance \U2 V — C/f°| as a function of j y . A least square fit on the last 
computed values (in log-log scale) gives \U^ y — £/f°| ~ Jy 2 ' 6 - 

We present in Figure 14.61 the dependence of the viscosity on the friction param- 
eter 7 y , for a fixed value j x = 1. A mild dependence on ^ y is observed in the limit 
j y — > 0. The value obtained for j y +oo is on the other hand very different from 
the limit obtained as j y — > 0. Note also that the viscosity seems to be an increasing 
function of the transverse friction, which makes sense from a physical viewpoint. 

4.2.4. Dependence and asymptotics in j x . The discussion after Theorem[3] 
suggests that the longitudinal velocity decreases as l/j x as j x — > +oo. To observe 
numerically this behavior, it is necessary to increase the magnitude of the nongradient 
force. Otherwise, the response is very small (indeed, proportional to r ) x 1 ) and relative 
statistical errors are too large to obtain meaningful results. We therefore computed 
the linear response of the velocity for values of £ proportional to -f x . This is done by 
modifying the evolution on p x in (|2.1j) as follows: 



d Pxi ,t = -V tol 7(ft) dt - lx _ CF(q yi , t )) dt +Jjf dW t x \ 
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Fig. 4.5. Convergence of the velocity profile for increasing values of the transverse friction ■jy. 
The dashed line represents an affine fit in log-log scale. 
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Fig. 4.6. Shear viscosity n as function of ^ y in the case ■yx = 1, for the sinusoidal nongradient 



This amounts to replacing £ in (|2.ip by £ = 7^^. The resulting average velocity profile 

IS 7a; ^a; ■ 

The results depicted in Figure l4~7l show that the average velocity, properly rescaled 
by j x , converges to the nongradient force, as predicted by ()3.15j) . The estimated 
convergence rate is I'YxUi" — ~ lx °' 9 - 

The behavior of the corresponding viscosities cannot be predicted from the results 
of Theorem[3] (see the discussion at the end of Section l3.4.2p . We therefore investigated 
numerically this dependence, see Figure 14.81 The viscosity is more or less constant 
for low values of j x , and increases for larger ones. 



16 



R. Joubaud and G. Stoltz 



<XX3K 



1 , 

8 

b 



















~~x_ 






~ x -x. x 















7* 



Fig. 4.7. Convergence of the velocity profile for increasing values of the friction -y x . The dashed 
line on the right picture represents an affine fit in log -log scale. 




Fig. 4.8. Shear viscosity n as function 0/7^ in the case -f y = X, for the sinusoidal nongradient 
force. Left: behavior for small values of "f x . Right: large 7^ asymptotics. 



5. Proof of the results. Unless otherwise stated, the norm || ■ || refers to the 
norm induced by the canonical scalar product on L 2 (t[)o). The operator A a .thm (fit — 
x,y), defined in (I3.10[) . can be rewritten as 

1 N 

i=l 

Note also that 

m 

where [A, B] = AB — BA is the commutator of two operators. 

5.1. Proof of Theorem [TJ The existence and the uniqueness of the invariant 
measure which has a smooth density with respect to the Lebesgue measure for any £ S 
M is a standard result since the position space is compact and the forces are smooth. 
It suffices to use hypoellipticity arguments and take the kinetic energy as a Lyapunov 
function (see for instance [35] for the general strategy, and [27J Appendix A] for the 
specific case under consideration). As a consequence, and recalling the definitions of 
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the operators given in Section 13.11 

Ker(Lo) = Span(l) = {cipo, c G K} , 

the vector space of constant functions on L 2 (-0o). Note also that Ker (Aq) = Span(l) 
by Proposition 15]. 

The key result to prove the expansion (I2.3[) is the following lemma (proved below). 

Lemma 1. The operators Lq 1 and Lq 1 B* are bounded operator on H (endowed 
with the L 2 (ipo) scalar product). 

In view of this result, we can introduce 

and define, for k > 1, 

ffe+i = — L 1 £Tffc, 

with 

h = -Lq 1 (6*1) = -lh 1 fepeiFidyi)) , (5.1) 

m \i=l / 

which is well-defined since the function (q,p) H> p X iF{q y i) belongs to % for all i = 
1, . . ,,N. The function in f|2.3f) is well defined for |£| < £* and a straightforward 
computation shows that 

Ufr = o- 

The uniqueness of the invariant measure allows to conclude. 
We now write the 

Proof of Lemma [7J We denote by || • || the L 2 (^>o)-norm. Standard results of 
hypocoercivity show that Aq 1 is bounded on W (and in fact compact, see [32l 14} [T3l 
[25] ) . Besides, for a smooth test function ip, 

N 

\\B V \\<\\F\\ L ~J2W d p*Ml 

i=l 

while 

(<p,Ao<p) = -~ (7x||V p ^|| 2 +7 S ||V P ^|| 2 ) . 

This shows that there exists a constant C > such that, for any smooth test function 

f, 

\\Btp\\ 2 <C\(<p,Ao<p)\<C\\<p\\\\A 'p\\. 
In conclusion, for any tp G H, 

IIBA-Vll 2 < c\\a^M |M|. 

Since Ran(£>) C "H, this shows that is a bounded operator on %. The same 

holds true for £ ~ 1 B*\-h, which is its adjoint on % C L 2 (ipo). □ 

Remark 1. In the above proof the fact that both j x and 7 y are non-zero is a 
crucial assumption. If for instance j x = 0, many arguments break down, and the 
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proofs become much more technical and/or some results cannot be proved anymore. 
This is due to the fact that the the Lie algebra generated by {Ahunn d Pyl d PyN } may 
be different from the Lie algebra generated by {Aham, d Pvel , • • • , 9 PxN , d Pyl , . . . , d PyN }. 
Indeed, 

N N 

[A am , d p J = d qyt , [A am , d q J = J2 C«* v ■ d Pvj + J2 v ■ d p*i ■ 

Possibly, iterated commutators should be computed as well. Additional assumptions 
on the potential are required to infer that d Pa:i is in the Lie algebra. This amounts to 
assuming that the coupling between the x and the y directions is strong enough. To 
our knowledge, the only cases where such arguments could be used are one- dimensional 
atom chains, for which the simple geometric structure of the system is of paramount 
importance to show that the Lie algebra has full rank. Obtaining hypocoercivity esti- 
mates is more challenging and imposes further restrictions on the interactions (see JS/). 



5.2. Proof of Proposition [TJ Corollary Q] shows that 
km 



r Wggjj>C P / TTe(v „, 

J lm n c = \ U x{Y,q,p),2_^p xi F(q vi ) 



1=1 I £ 2 M>o) 



L N f 



where ipo(q) dq = Z~ x e^^ v ^ q ' dq is the marginal of the canonical measure in the 
q variable. The integrand in the last equation depends only on one variable q y i. By 
translation invariance of the system, the marginal distribution in the q y i variable of 
ipo(q) dq is the uniform distribution on L y T. Therefore, for any i = 1, . . . , N, 



Xe(q yi - Y)F{q yl ) i> {q) dq = f - / " X e(v - Y)F{y) dy — > F(Y) 

Lin 



1 

IV N Ly JQ 

as e — > 0, so that 



Um lim x l = --FCY). (5.2) 



Now, a simple computation shows that 

AaU*(Y, Q,P) = -^~ x (j2 ^^nX's (fc, -Y)- X * (?* - Y) d^V^ 
-^U s x (Y,q,p). 

The sum on the right-hand side can be decomposed into two contributions, one pro- 
portional to the kinetic part of the off-diagonal part of the stress tensor, and the other 
one arising solely from interaction forces. The first contribution can be written as 



d^ xvMn (Y,q,p) 



d ( 1 ^PxiPy t , v \ 



dY dY \pL x 4^ m 
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For the second part, we first use the pairwise character of the interactions to write 



d^V{q) = Y,V'{\q i -q j \) 



Qxi Qxj 

h -qj\ ' 

and then symmetrize the resulting expression as 

N 

£ Xe (% ( ~ Y) d q ^V(q) =Y,Xe (q vt - Y) Vflft - 9i\)^H, 
i ■ i - 

= E (Xe(qyi-Y)-Xe(qy j -Y))V'(\q i -q j \)^^f. 

l<i<j<N * 

The second contribution finally reads 
dZ E pot (Y,q,p) d I 1 ^ ^fQxi-QxA f qvt , v ,. 



In conclusion, it holds 



AoU £ x (Y,q,p) = "" v ~ —UZ<y,q,p). (5.3) 

Combining the latter result and (I5.2[) leads to 

1- li- i f ^ (S y 0>g + ^ (^(y, 0>^ = -FiX). (5-4) 



■/?; 



Now, Corollary Q] shows that the limit 

7 ;!V--).- . i„ x 

»=i 



<4P0 : = J im — ^-1— = V / PxiXs(q V i -Y){i(q,p)tpo(q,p)dqdp 



is well defined. By hypoellipticity, the function fi belongs to C°°(V N x R 2Ar ). The 
limit e — >■ of the right-hand side is therefore well defined and 

L N f 

u x (Y) = -t^-Y] / p xi (f l ip )(qi, . . . ,qi-i,q xi ,Y,q i+1} . . . ,q N ,p) dqi-.i-.N dq xi dp, 

i=l Jv N - 1 xL :c TxR 2N 

where dqi.i-.N = dqi . . . dqi-i dqi+i . . . dqN- A similar reasoning holds for <r xy . Passing 
to the limit in (|5.4p . 

p oY m m 

which is (13.61). 



5.3. Proof of Theorem [2j To simplify the notation, we set m = 1 in this 
section, but the proof can be straightforwardly modified to account for more general 
masses. Note first that the solution of p. lip is well defined for any j y > by the 
Fredholm alternative (since -4o(7j/) nas a compact resolvent on H = L 2 (ipo) n {1} , 
and the right-hand side of the equation is orthogonal to Vect(l)). 

We start by formal computations providing possible expressions of / , / , and 
then prove rigorously the convergence result stated in Theorem [5J To this end, we 
need some intermediate uniform hypocoercivity result. 
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Formal asymptotic expansion in "f y . We consider the following ansatz for the 
solution / 7 : 



lv % 



/7„=/° + -/ 1 + 4/ 2 



and rewrite the operator Aq{^ v ) as the sum -4o(7y) = To + J y Ay, thm- The kernel of 
the operator A y ^m on L 2 (-0 O ) is 

Ker(A Vt tbm) = {.9 6 £ 2 (V>o) 5 = s(?,Px)}- 
This is a consequence of the equality 

(g,A y ,thmg) L 2 {M = -^I|V Ph 3|| 2 

and the fact that the Gaussian measure (in the p y variable) satisfies a Poincare inequal- 
ity. Identifying terms with the same powers of *y y in (|3.11l) , the following hierarchy is 
obtained: 



N 



T f° + Ay^hmf 1 = -^PxiG(q y i), 

i=l 

> Tof + Ay^thmf 2 = 0. 



(5.5) 



The first equation shows that f° = f°(q,p x ). The second equation can then be 
rewritten as 



JV 



Ay^ m f X {q,p) = -p y ■ ^7 qy f°(q,Px) - ^2pxiG(q yl ) - T qy f°(q,p x ). 



where 



T qy =p x ■ - V qx V(q x ,q v ) ■ V Pm + j x Aa 



,thn 



is an operator parameterized by q y G (L y T) N , and acting on the Hilbert space 
L 2 (^ qv ), where 

% v (q x ,p x ) = Z^exp (-0 (v(q x ,q y ) + f^)) ■ 

Setting 

f 1 =J 1 +PyV qv f, 



it holds 



A y ,thmf 1 (q,P) = -^P X iG{q yi ) -TqJ°(q,p x ). 

i=l 

Since the right-hand side does not depend on p y , the solvability condition for this 
equation is that the right-hand side vanishes. Besides, by the same results which allow 
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to prove that Aq has compact resolvent, the operator 7~ q , considered as an operator 
on L 2 (^ q ), has a compact resolvent on Ker(T qy ) ± = {1}^ (where the orthogonality 
is with respect to the canonical scalar product on L 2 {^ q ); see (351 Proposition 15] 
for a proof of the latter equality). In addition, by hypocoercivity, T qy x is bounded on 

H 1 (9q y ) nKer(7^ y ) ■ Therefore, T q ~ 1 {p X i) is well defined by the Fredholm alternative. 
By linearity, 



N 



f(q,p) = -Y J G{q yl )T qy 1 (p xl ), 



i=l 



and 



f 1 (q,p)=PyV qy f(q,p) + f 1 (q, Px ), 



provide admissible solutions for the first two levels of the hierarchy (|5 . 5|) . The func- 
tion f 1 will be made precise below (see (|5.12p V 

The function f° is in H 1 ^®). To show that the function f 1 — f 1 is indeed well 



defined, it is enough to show that d qyi T q (p x k) 



is well defined. This, in turn, 



follows from the following equality for any function ip = Lp(q x ,p x ): 



N 



if. 



(5.6) 



The operators d Pxj [Tq^j are bounded on L 2 (^ qy ) n {1} for j = 1, . 

(rt^W,,) = -j\\VpM\l*(* qy) , 

which implies, for ||v 5 IIl 2 (*, ) _■ 1; 



, N since 



•r 



L 2 (^ qy ) l x 



< 



(v) 



L 2 (-f qy ) l x 



'In 



In conclusion, J 1 — J 1 € H 1 ^^). In addition, by hypoellipticity, the functions /°, J 1 — 
J 1 are C°° when G is smooth. 

Uniform hypocoercivity estimates. Let us show that the operator Ao("j y ) is uni- 
formly hypocoercive for j y large enough (say, j y > 7 X ), provided the domain of the 
operator is restricted to functions with vanishing average with respect to the Gibbs 
measure in the p y variable. To this end, we decompose Aoi'Jy) as 

A>(ly) = Milx) + {l y ~ Jx)Ay,th m . 

Following the proof of Theorem 6.2 in [11], it can be shown that there exists k > 
such that, for all smooth functions u € %, 

- ((u,Aq(j x )u)) > k((u,u)) , 

where the norm induced by ((-,■)) is equivalent to the i? 1 ('0o) norm 



H 5 



IVoUll 



|V,u| 
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More precisely, ((■,■)) is the bilinear form defined by 

{{u, v}} = a{u,v) + b (V p u, V p w) + (V p u, V q v) + (V 9 u, V p u) + b(V q u, V g w), 

with appropriate coefficients a » & s> 1 ■ It follows that there exists C > indepen- 
dent of 7 y such that 

C'HIiUwo) ~ (7y -7a:) ««,A/.thm«)) < ~ -^o(7» )«)) ■ 

Let us now show that 

- {{u,Ay, thm u)) > (5.7) 

for functions u in an appropriate subspace of i? 1 ( - 0o)- Using the commutation rela- 
tions [d Pa , .] = f36 ata >5ij (a, a' € {a;, y}), a simple computation shows 

=£(« + ^)ll^ v ^ll 2 + &llv p a Pv ^|| 2 

+ &|| V q d Pyl u\\ 2 + 2(V q d Pyl u, V p d Pyl u) + P(d qyi u, d Pyi u) 

N 



> E ( a + ( b - \ ) ) h^ii 2 + ( 6 - v ^«n s 



+ (&-i)l|v 9 a PyiU || 2 -^||5 gHiW || 2 . 

Summing on i £ {1, . . . , N}, the quantity (|5.7I) is seen to be non- negative for an 
appropriate choice of constants a 3> b ;§> 1 provided there exists a constant A > 
such that, for alH = 1, . . . , N, 

\\d qyi u\\ <A\\V p d qvi u\\- (5-8) 

This indeed implies 

N N N N 

E K^u 2 < A E ii^a^ii 2 = ^E n v A ii 2 < ^E i! v A, ii 2 - 

i=l i,j=l j=l j=l 



Since the Gaussian measure satisfies a Poincare inequality, the inequalities (|5.8[) hold 
provided 

Vi = 1, . . . ,N, / d qyt u(q,p) cxp I -j8-^ ] dp y = 0. 



Defining the closed subspace of L 2 (ipo) H {!}" 



H = «eflU) 



«(?,Px) = ( -J 



2^ - W/2 



jf ^ v(q,p) exp (~0yj - I C H, 



(5.9) 

we conclude that, for any ueWofl H 2 (tpo), 

C\\u\\ 2 HHM <-({u,A { ly )u)). (5.10) 
In particular, there exists a constant K > such that, for any 7 y > 7^ and for any 

||A(Ty) _1 w|| fflw , o) < K\\u\\ H i^ o) . 
In fact, this inequality can be extended to functions in %q. 
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Proof of the limit (|3.12l) . To prove (|3.12p . we proceed as follows. Note first that 

-Mlv) {fly -f~ = -Tof 1 , 

so that 

ffy ~ f ~ ly'f 1 = --M^Tof. (5.H) 

Since „4o(7j/) _1 is bounded on T-Lq, uniformly in j y (see (|5.10p ). it is sufficient to show 
that Tof 1 6 Ho- In view of the definition of /°, the proof is then concluded by setting 

<t>i(l'P) = -Tq^iPxi)- 

Let us first show that Tof 1 (q,p x ) = (where v is defined in (JOS))- This can be 
ensured by an appropriate choice of / . Note first that 

T f = P y • VgJ 1 + TgJ 1 + ( P y • V q y - V q y V ■ V Py ) f 1 + T q J\ 

The first two terms have a vanishing average with respect to (2*)"^ exp dp 
Introducing 

g(q,p x ) = -(2n)- N / 2 f ( Py ■ \7 qy - V qy V ■ V P J f 1 exp ( J d Py , 



v 



K N 



the condition Tof 1 = is satisfied provided 

Seeing the function on the right-hand side as a function of (q x ,p x ) indexed by q y 
allows to define f 1 pointwise in q y as 

f 1 = -{2*)- N l*T- l g. (5.12) 

Let us now study the regularity of To/ 1 . We only treat the term Toif 1 — f 1 ) 
since the regularity of To/ 1 can be proved similarly. Recall that all the functions 
under consideration are C°° by hypocllipticity. Therefore, only the derivates in the p 
variables have to be considered because the position space is compact. Now, 

f 1 -7 1 =-^PyiG'(Q«i)7^ 1 M 

i (5 13) 

- ^2p yi G(q yj ) {(T- 1 ) [d 2 qv ^ k V(q x ,q y )d Pxk ] (V^ 1 ) Pxj ) . 

The p y dependence is trivial in the above expression, so that only derivatives in p x re- 
quire some attention. Since T = A Vt \ xa , ra +T qv where -4 y ,ham = Py • V 9y - V qy V(q x , q y ) ■ 
V Pb is an operator in the q y ,p y variables (parameterized by q x ), it suffices to consider 
T q f 1 . This function is, in turn, a linear combination of terms of the form p y ip X i (cf. 
the first term in the right-hand side of (|5 . 13|) ^) and Pyid P:ck T q ~ 1 p X j (second term in 
the right-hand side of (|5.13p V To prove that the latter functions are in i? 1 ('0o), we 
use the results of [32JEI], which show that T^ 1 is a bounded operator on the Hilbert 
spaces 

\feH m (V qy ) [ f(q x ,p x )* qy (q x ,p x )dq x dp x 1 Ci 2 (\) 

[ J(L m T) N xR N J 

for any m > 0, with a bound uniform in q y . 
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5.4. Proof of Theorem [3j The proof follows the same lines as the proof pre- 
sented in Section 15.31 so we skip the parts of the argument which can be straightfor- 
wardly extended from there. 

As in the previous section, we set m = 1 to simplify the notation, but the proof 
can be straightforwardly modified to account for more general masses. Note first that 
the solution of (|3.13j) is well defined for any -f x > 0, for reasons similar to the ones 
exposed at the beginning of Section 15.31 Define 

Tq x = p y ■ V qy - V qy V(q x ,q y ) ■ V ' Py + J y A y , t hm, 

which is an operator parameterized by q x S (L X T) N , and acting on the Hilbert space 
L 2 (^ q J, where 

y qx (q y ,p y ) = Z->xp ^-/3 (v(q x ,q y ) + %L Yj . 

Its kernel is Vect(l) = {c* gx , c e M}. 

Formal asymptotic expansion in j x . We start by formal computations, with a 
discussion parallel to the corresponding one in Section 15731 We consider the following 
ansatz for the solution f~ a : 

/ 7 ,=/° + -/ 1 + ^/ 2 + ... 

and rewrite the operator Ao(~f x ) as the sum Ao(~f x ) = Tq + j x A Xt thm- The kernel of 
the operator Ax,tbm on L 2 (ipo) is 

Ker(A,thm) = {9 e ^ 2 (V'o) 9 = g(q,p y )\- 

Identifying terms with the same powers of -f x in (|3.13[) . the following hierarchy is 
obtained: 

A,thra/ 0> 

N 

T f° + A.thm/ 1 = -^PxiG(^), (5.14) 

i=l 

, Tof 1 + A x ,thmf = 0. 

The first equation shows that f° = f°(q,p y ). The second one can be rewritten as 



so that 



N 

A^thm/ 1 (<?,£>) = -^2pxt(G(q yi ) +d qxi f°(q,py)) - T qm f°(q,p v ), 



N 

/ 1 = &**( G ^)+^*/V^))+/ 1 > 



with 



Ac.thm/ 1 ^^) = -Tq x f°{q,Py). 
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The solvability condition requires T qn .f°(q,p y ) = 0, hence f° = f°(q x ). Besides, f 1 
does not depend on p x . The solvability condition for the third equation in (|5.14[) is 
To/ 1 G Ker(A,thm). Now, 

(N \ N 

X>-(G(9 V i) + ^/°(g x )) =-^9 toi n&,%)(G(^)+^ i /°(&) 
i=l / i=l 

N N 



1=1 



and 



i=l 



AT 



We therefore set 

/° = 0, f 1 = T,- 1 9, xi V(q x ,q y )G(q yi )j , 

and 

N 

f 2 (q,p) = ^p„(p st G'(? 51 ) + a^j 1 ) + f 2 , 
1=1 

so that 

TV 

A,thm/ 2 (g,p) = - ^Pxi^tfiG'^i) + dq^f 1 ) = -To/ 1 . 

i=l 

The function f 2 is chosen such that To/ 2 has a vanishing average with respect to the 
Gaussian measure in the p x variable. 

Note that f 1 is well defined since Yli=i 9q xi V(q x , q y )G(q y i) G Ker(T 9x ) ± = {l} ± 
(where the orthogonality is with respect to the scalar product on L 2 (\I> qx )) . Indeed, 

N N 

^d qxi V(q x ,q y )G(q yi ) = J2^2V(\q i -q j \)^^iG(q yi ) 

i=i i=ijy« iqi qjl 

= £ V'(\ qi - qj \)^^f(G(q yi )-G(q yj )). 

l<i<j<N W ®l V ' 

The last line is antisymmetric with respect to the exchange of coordinates q y t and 
q y j , hence the average of the corresponding function with respect to 4^ , which is 
symmetric with respect to the exchange of coordinates q y i and q v j, vanishes. 

A discussion similar to the one in Section [5731 show also that c^/ 1 is well defined, 
hence the definition of f 2 makes sense. 

Proof of the limit (|3.14[) . The remainder of the proof follows the very same lines 
as the proof presented in Section [5. 31 hence we omit it. 
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